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Resume. Nous proposons une nouvelle methode pour estimer les ruptures du rythme cardiaque dans les 
bandes parasympathiques et orthosympathiques basee sur la transformee en ondelettes dans le domaine 
complexe et P etude des ruptures des moments des modules de ces transformees en ondelettes. Nous obser- 
vons des ruptures dans la distribution pour les deux bandes de frequence. 



1 Introduction 

ECG signal analysis has a long story after the implementation of the ambulatory monitoring by Holter at the beginning 
of the fifties. Recent measurement methods, due to the size reduction of the measurement devices, see Chamoux (1984), 
allow us to record ECG data for healthy people over a long period of time : long distance (marathon) runners, individuals 
daily (24 hours) records, etc. We then obtain large datasets that allow us to characterize the variations of the heartbeat rate 
in the two components of the nervous autonomous system : the parasympathetic and the orthosympathetic ones. 




FlG. 1 - RR interval for a healthy subject during a period of 24 hours 
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The data studied in this paper have been recorded the 29 of January 2008, at the Clermont-Ferrand University 
Hospital (CHU) by G. Boudet. Unlike Boudet et al. (2004) and Cottin et al (2006), these data are not recorded in the 
framework of a laboratory experiment, but during the real life of the individual under investigation. 

We have daily datasets, over 100,000 observations, and in the near future we will have data over several days, i.e., over 
that sample size. We are measuring the interval between two RR peaks, i.e., if we denote by (tj)j=i,...,iv, the sequence of 
peak times measured with a precision of 1 .E-03 second, we consider the series X (t , ) = (tj — i ^ measured in seconds. 
Several indicators are related to heartbeat data. The most popular is the instantaneous average frequency, i.e., Xft.i)^ 1 , 
which is displayed by runners watches. This quantity is informative on daily activity of observed individual : sleeping and 
waking up times, physical activity, e.g., sport, manual work, but does not summarize all relevant information. Note that 
heartbeat data display large variations, clustering, etc, as only individuals with serious disease display a regular heart rate. 

Cardiologists are interested in the study of this signal in two frequency bands : the orthosympathetic and parasympa- 
thetic bands, i.e., the frequency bands (0.04 Hz, 0.15 Hz) and (0.15 H z, 0.5 Hz) respectively. The definition of these 
bands is the outcome of research works, see e.g., Task force of the European Society of Cardiology and the North Ameri- 
can Society of Pacing and Electrophysiology (1996), and is based on the fact that the energy contained inside these bands 
would be a relevant indicator on the level of the stress of an individual. 

Indeed, for the heart rate, the parasympathetic system is often compared to the brake while the orthosympathetic 
system would be a nice accelerator; see e.g. Goldberger (2001). At rest there is a permanent braking effect on the heart 
rate. Any solicitation of the cardiovascular system, any activity initially produces a reduction of parasympathetic brake 
followed by a gradual involvement of the sympathetic system. These mechanisms are very interesting to watch in many 
diseases including heart failure, but also rhythm disorders that may fall under one or other of these two effects, monitoring 
the therapeutic effect of several Medicines including some psychotropic. In the field of physiology such data are crucial 
for measuring the level of vigilance and particularly the risk of falling asleep while driving a vehicle, the level of stress 
induced by physical activity or level of perceived stress, which can be considered as a criterion of overtraining in sport. 

Fractals models have been used in cardiology after the works by Ivanov et al. (1999), who applied the multifractal 
spectrum analysis advocated by Frisch (1996), for modeling RR series and classifying individuals according to this mul- 
tifractal spectrum, as this spectrum discriminates between individuals who experienced hearth trouble, and those who did 
not. However, this tool has some shortcomings as it requires huge samples common in turbulence analysis, and is then 
inappropriate for studying phenomena occurring at a resolution lower than the daily time interval, such as the variations 
of the parasympathetic and orthosympathetic systems inside the day. 

Wavelets based methods have been used in biostatistics by Diab et al. (2007) for uterine EMG signal analysis. Ho- 
wever, they consider that the process studied is homogeneous, and used these methods as a classification tool. Another 
significant difference is the fact that they use discrete wavelet decomposition, i.e., a frequency decomposition on a dyadic 
wavelet basis, the choice for the frequency bands is made without reference to a biological phenomenon. In our case, the 
choice for the frequency band is justified by biological considerations, and we fit the wavelets inside these bands. 

This is why the continuous analysis of both systems and their quantification is a particularly promising research area. 
So, the example in Figures 4, 5, 6 shows at the observation 28,220 i.e., 13h40'50", a simultaneous variation of both 
systems. 

2 Mathematical formulation 

For biological reasons, the heart rate is within the interval [20, 250] bpm (beats per minute), i.e., X(t) belongs to the 
RR interval 60/250s. < X(t) < 60/20s.. This leads to modeling by stationary or locally stationary processes. We wish 
to decompose this signal in a sequence of homogeneous intervals having the same mean and/or variance. We make the 
assumption that the series is Gaussian, and then find the optimal segmentation of the process using the same approach as 
in Lavielle and Teyssiere (2006) or Lavielle and Moulines (2000). 
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We assume that the signal is the sum of a piecewise constant function and a Gaussian process, centered and locally 
stationary. We then have the following representation : 

X{t)=p(t)+ [ ^f 1/2 {t,^)dW{^), forall t e R, (1) 

J Ft. 

where 

- ^ i ► f(t, £) is an even and positive function, called spectral density piecewise constant, i.e., there exists a partition 
Ti, . . . , t k such that f(t, £) = f k {£_) for t 6 [7*, n +1 [ 

- the function t 1— > p(t) is also piecewise constant for another partition tx, ■ ■ ■ , Tr, with p(t) = /i£ if t £ [t£, Tf +1 [ 
The wavelet coefficient associated with ip is 

W^{b)= I tp(t-b)X(t)dt unit in second 2 . (2) 

In Bardet and Bertrand (2007) or Bardet et al. (2008), one can find a theoretical study of the wavelet coefficient for 
stationary (or with stationary increment) centered Gaussian processes, i.e., for X given by ([]]) with p(t) = 0. This 
approach can be generalized to locally stationary Gaussian processes, this will be detailed in a forthcoming paper. 

According to recommendations of Task force of the European Soc. Cardiology and the North American Soc. of Pacing 
and Electrophysiology (1996), we use the following notations : 

- [wi,a->2] = (0.04 Hz, 0.15 Hz) denotes the orthosympathetic frequency band, 

- [u>2, CJ3] = (0.15 Hz, 0.5 Hz) denotes the parasympathetic band 

The energy associated with each of these frequency bands and localized around the time b, is measured by the modulus of 
the complex wavelet coefficients |Wj(&)| 2 for i — 1,2, with 



b\ — >Wi{b) = / 4>i(t-b)X(t)dt, i = l,2. 

JFt 

In this work, these wavelets coefficients are computed at each second, i.e., the difference between two consecutive values 
for b is equal to 1 second. 

How to choose the wavelets ipi and 1P2 ? 

In the idealistic case, one would use two filters ipi and 1P2 with compact support, the Fourier transforms of which have 
support inside the orthosympathetic and parasympathetic bands. 

Unfortunately, it does not exist a non null function tp with compact time domain support and compact frequency 
support, see for instance Mallat (1998) Th 2.6 p. 34. Therefore, the best we can do is to choose between a filter with a 
compact frequency support and a filter with a compact time domain support. The first choice is well suited for stationary 
models, see Bardet and Bertrand (2007). But, in this work, we are interested by locally stationary models, thus our 
specifications are a filter with compact time domain support as a Daubechies wavelet. The price to pay for the compactness 
of the time domain support is the loss in the compactness of the frequency support. However, the frequency support is 
"almost compact" in the following sense : 

Definition 2.1 (p pseudo support) Let < p < 1, be a map g E L 2 (M) that admits the compact interval I as a p pseudo 
support if = p. 

J mi \g{t)\ 2 dt 

Fourier transform of Daubechies wavelet have a reasonably small p pseudo support with a ratio p close to 1. Moreover, the 
larger the number of the Daubechies wavelet is, the closer to 1 the ratio p is ; see the example below with the Daubechies 
wavelet D6. 

By scaling and modulation, one can adjust the pseudo support inside a specified frequency band as stated by the 
following proposition 
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Proposition 2.1 Let ip be a filter with compact support [L\, L 2 ] and a frequency p pseudo support [Ai, A2], for any 
frequencies band [coi, UJ2], the map ipi(t) = fj, X e lr,t ip(\t) with 

W2-W1 , LU1+UJ2 , . A 2 + Ai 

p > 0, A = — and t] = (lu 2 — ui)- 



A 2 - Ai 

has a p pseudo support [uj\ , a-> 2 ] and a time domain support 



A 2 - Ai 



A 2 - Ai A 2 - Ai 
M, 1^1 

U)2 — UJ\ W 2 — U>\ 



Proof. Since V'i(£) = A* x "0 
proposition. 



A 



one can deduce p pseudo supp ipi = 77 + A x p pseudo supp ip and then the 



Different choices for the wavelets ipi and V'2 



From Proposition ^. 11 one can deduce the different possible choices of the filters ipi and ip 2 

- Daubechies wavelet D6 : In this case, we have A x = 0.08, A 2 = 1.75, p = 0.9999 and one can set 

ip 1 (t) — fix e imt D e (\it) and tp 2 (t) = p x e 1 ^ 1 D 6 (X 2 t) 
with ??! = -0.0255, Ai = 0.0659, 7 ?2 = -0.0585, and A 2 = 0.2096. 

The length of the time support are \Suppipi\ and \Supptp2\- One can see on Fig. 2 (a) that the Fourier transforms 
•tpi{x) and 4>2{x) have almost disjoint supports. 

u)\ + ll>2 , s A 2 + Ai 

Aj = — and Vl = (wa - ^i) A9 _ Ai 



0J2 


- UJl 


A 2 


-Aj 




- L02 


A 2 


— Ai 



u) 2 + u 3 A 2 + Ai 

A 2 = — and i]2 = g ^ 3 ~ ^) Az _ Ai 

\Supptpi\ = — x \SuppDe\ and \Supptp2\ = —— x \SuppD§\ 

0J2 — u)\ 0J3 — ui2 

Gabor wavelet : For computational reasons, we will use the Gabor wavelet defined as 

m = e^9(t), 9 (t) = J -± IFi e-£ (3) 

see, e.g., Mallat (1998). This wavelet has the same time and frequency p pseudo support [-L, L] = [—3.5, 3.5] with 
p = 0.9995. In the spectral domain, we have 

m=9(Z-v), g(0 = (^ 2 ) 1/i e^ (4) 



We fit the Gabor wavelet inside the orthosympathetic or the parasympathetic frequency bands by using Prop. 12.1 l or 
direct calculations. One obtains the two Gabor wavelets defined by (0 with the following choice for the parameters : 

u)\ + LO2 2L 
7/1 = and 01 



2 LO2 — u>\ 

L02 + W3 2L 
Tj2 = ^ and fj 2 



2 LO3 — UJ2 

4L 2 4i 2 
moreover \p pseudo Supp = and \p pseudo Supp il> 2 \ = with p = 0.9995 

One can see on Fig. 2 that the Fourier transforms (x) and i/) 2 (x) still have almost disjoint supports. 
Using the Gabor wavelet is more efficient in terms of computing time, as it is at least 8 times faster. Figure 3 below 
displays the Gabor wavelets coefficients in the orthosympathetic and parasympathetic bands respectively for the sample 
plotted in Figure 1. 
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FlG. 2 - The Fourier transforms ipi(x) (left) and ip2{x) (right) 
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FlG. 3 - The wavelet coefficients in the orthosympathetic band (left) and in the parasympathetic band ( right) of the same 
healthy subject during a period of 24 hours 



3 Segmentation analysis 

We assume that the process {X t } is abruptly changing and is characterized by a parameter 9 € that remains constant 
between two changes. We use this assumption to define our contrast function J(r, X). Let K be some integer and let 
t = {ti, T2, . . . , tk-x] be an ordered sequence of integers satisfying < t\ < < . . . < tk-i < n. 



For the detection of changes in the mean and variance of a sequence of random variables, i.e., a change in distribution, 
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we use the following contrast function, based on a Gaussian log-likelihood function : 

J n (r, X) = - ]T Ut " :* tJ2 + n k log(^). (5) 

where n k = — Tk-i is the length of segment k, a\ is the empirical variance computed on that segment k, <j\ = 
n k 1 SI=r fc 1+1 (-^i — X) 2 , an d i s me empirical mean of X\, . . . , X n . 

In fact, we minimize a penalized version of this contrast function, the penalty parameter being selected in an adaptive 
way so that the obtained segmentation does not depend too much on the penalty parameter ; see also Birge and Massart 
(2007) for another choice for the penalty parameter. Although this method has been devised for Gaussian processes, it 
empirically provides relevant results for non-Gaussian processes, e.g., financial time series, see Lavielle and Teyssiere 
(2006) for further details. 

3.1 Orthosympathetic band 

For that frequency band, we obtain the following segmentation : 

r = {28220,33366,71048}, 

i.e., 13h40'50", 15h06'36", lh34'38". 




FlG. 4 - Segmentation in the mean and variance (change in the distribution) of the modulus of the wavelet coefficients in 
the orthosympathetic band 

3.2 Parasympathetic band 

For that frequency band, we obtain the following segmentation : 

t = {11620, 21912, 28054, 31540, 36022, 40172, 52622, 70550}, 
i.e., 9h04'10", llh45'42", 13h38'04", 14h36'10", 15h50'22", 17h00'02", 20h27'32", lh26'40". 

4 Conclusion 

The example illustrated by Figures 4, 5 and 6 shows at around t = 28, 220 a simultaneous variation of both systems. 
But, one can also observe change-points in the orthosympathetic and parasympathetic bands occurring at different times. 
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FlG. 5 - Segmentation in the mean and variance (change in the distribution) of the modulus of the wavelet coefficients in 
the parasympathetic band 



In the future, we will study the existence of a possible causality or sequentiality between these change-points in different 
bands. 



FlG. 6 - Mean of the modulus if the wavelets coefficients for the orthosympathetic and parasympathetic bands, and of the 
RR process X(t) 
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Summary 

We propose a new method for estimating the change-points of heart rate in the orthosympathetic and parasympathetic 
bands, based on the wavelet transform in the complex domain and the study of the change-points in the moments of the 
modulus of these wavelet transforms. We observe change-points in the distribution for both bands. 



